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ABSTRACT 

A  fast  algorithm  for  recovering  profiles  of  density  and  Lame 
parameters  as  functions  of  depth  for  the  inverse  seismic  problem  in  an 
elastic  medium  is  obtained.  The  medium  is  probed  with  planar  impulsive 
P  and  SV  waves  at  oblique  incidence,  and  the  medium  velocity  components 
are  measured  at  the  surface.  The  interconversion  of  P  and  SV  waves  defines 
reflection  coefficients  from  which  the  medium  parameter  profiles  are 
obtained  recursively.  The  algorithm  works  on  a  layer-stripping  principle, 
and  is  specified  in  both  differential  and  recursive  forms.  A  physical 
interpretation  of  this  procedure  is  given  in  terms  of  a  lattice  filter, 
where  the  first  reflections  of  the  downgoing  waves  in  each  layer  yield 
the  various  reflection  coefficients  for  that  layer.  A  computer  run  of 
the  algorithm  on  the  synthetic  impulsive  plane  wave  responses  of  a  twenty- 
layer  medium  shows  that  the  algorithm  works  satisfactorily. 
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INTRODUCTI ON 

A  layered  elastic  medium  with  depth-dependent  material  parameters 
is  probed  with  impulsive  plane  waves  at  oblique  incidence.  The  medium  is 
assumed  to  support  the  propagation  of  seismic  waves ,  so  that  there  is 
continuing  conversion  between  P  waves  and  SV  waves  as  the  medium,  with  its 
depth-varying  parameters,  is  penetrated.  The  goal  is  to  recover  profiles 
of  the  density  p(z)  and  Lame*  parameters  A  (z)  and  y(z)  as  functions  of 
depth . 

Previous  work  on  this  problem  has  yielded  methods  of  solution  that 
are  computationally  arduous  to  implement.  For  example,  Coen  (1983)  solved 
this  problem  by  employing  solutions  to  the  acoustic  problem  for  the  separate 
cases  of  P  and  SV  impulsive  plane  waves  at- normal  incidence,  which  are 
decoupled  for  a  layered  medium,  and  of  SH  waves  at  oblique  incidence. 

This  allowed  the  recovery  of  the  parameter  profiles  by  solving  Marchenko 
integral  equations,  but  sidesteps  the  issue  of  P-SV  mode  conversions. 
Blagoveschenskii  (.1967)  exhibited  several  integral  equations  whose  solutions 
yielded  the  parameter  profiles  as  functions  of  travel  times,  and  by  combining 
the  Gelfand-Levitan  inverse  scattering  method  with  the  solution  of  a 
Volterra  equation,  Carroll  and  Santosa  (1982)  were  able  to  recover  the  parameter 
profiles  as  functions  of  depth.  Baker  (.1982)  solved  the  related  problem  of 
reconstructing  radially  varying  parameters  by  using  spherical  harmonics  and 
Marchenko  integral  equations.  However,  none  of  these  methods  can  be  considered 
to  be  attractive  from  a  practical,  computational  perspective. 

Clarke  (1983)  and  Shiva  and  Mendel  (1983)  have  recently  given  algorithms 
that  utilize  the  layer-stripping  principle  employed  by  the  algorithm  given  in 
this  paper.  However,  their  algorithms  are  much  more  complicated  than  the 
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present  algorithm,  since  Clarke's  algorithm  requires  the  iterative  solution 
of  an  equation  at  each  step,  while  Shiva  and  Mendel's  algorithm  requires  the 
solution  of  a  cubic  equation  and  a  maximum-likelihood  estimation  at'  each  step. 
The  present  algorithm,  in  contrast,  is  quite  straightforward,  since  the 
assumption  of  a  continuous  medium  allows  differential  updates  of  the  medium 
parameters. 

In  this  paper  we  present  a  fast  algorithm  that  recursively  generates 
the  parameter  profiles  p (z)  ,  A(z)  and  ]i(z)  as  functions  of  depth.  There 
are  no  integral  equations  to  solve;  the  computation  is  direct  and  only 
involves  simple  operations.  The  algorithm  works  on  a  layer-stripping 
principle  that  involves  upgoing  and  downgoing  P  and  SV  waves  and  bears 
some  similarity  to  the  Schur  algorithm  and  the  downward  continuation 
method  which  were  used  to  solve  the  inverse  acoustic  problem  in  Yagle 
and  Levy  (1983)  .  and  Bube  and  Burri'dge  (1983)  ,  and  other  inverse  scattering 
problems  in  Dewilde,  Fokkema  and  Widya  (1981) ,  Bruckstein,  Levy  and 
Kailath  (1983).  and  Yagle  and  Levy  (1984)  ^  .  The  method  of  characteristics 
which  was  used  by  Symes  (1981) ,  Santosa  and  Schwetlick  (1982)  and  Sondhi 
and  Resnick  (1983)  to  reconstruct  the  impedance  of  an  acoustic  medium 
relies  also  on  a  similar  layer-stripping  technique. 

The  paper  is  organized  as  follows.  The  problem  is  set  up  in  detail 
in  the  next  section.  The  algorithm  is  derived  mathematically  and  exhibited 
in  the  following  two  sections.  The  next  section  discusses  what  the  algorithm 
is  doing  and  how  it  works  with  attention  given  to  physical  interpretations 
of  the  quantities  appearing  in  the  algorithm.  A  final  section  presents 
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the  results  of  a  computer  run  of  the  algorithm,  in  which  the  algorithm 
reconstructs  a  twenty-layer  medium  from  its  impulsive  plane  wave  responses. 

PROBLEM  FORMULATION 

We  consider  a  layered  elastic  medium  with  depth-dependent  density 

p(z)  and  Lame  parameters  AC z)  and  y(z).  The  following  experiment  is 

performed.  An  impulsive  planar  P  wave  is  obliquely  incident  on  the 

medium.  The  horizontal  and  vertical  velocities  are  measured  at  the 

surface  z=0  as  functions  of  time.  The  experiment  is  then  repeated  with 

impulsive  planar  SV  waves.  The  angles  of  incidence  of  the  P  and  SV 

plane  waves  with  respect  to  the  vertical  are  chosen  (i.e.  the  point  source 

data  is  stacked),  so  that  the  horizontal  ray  parameter  p  is  the  same  for 

both  experiments.  Of  course,  the  algorithm  may  be  run  concurrently  with 

many  different  values  of  p  from  a  single  point  source  experiment,  the 

updated  medium  parameters  at  each  depth  from  each  run  averaged,  and  the  averaged, 
values  then  used  in  the  algorithms.  This  reduces  the  effect  of  noise  in  the  data. 

There  is  a  choice  for  the  boundary  conditions  at  the  surface  (z=0) 

of  the  medium.  In  the  presentation  of  the  algorithm  a  free 

surface  is  assumed,  so  that  the  surface  tractions  T  and  T  are  both  zero. 

zx  zz 

Due  to  the  vast  difference  in  material  parameters  between  the  ground  and 
the  air,  this  assumption  is  quite  reasonable.  However,  the  medium  may  also  be 
considered  as  being  probed  from  an  overlying  homogeneous  half-space  in 
which  case  the  upgoing  P  and  SV  waves  at  the  surface  are  constructed  from  the 
velocities  at  the  surface.  Note  that  the  assumption  of  planar  source  waves 
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is  often  reasonable  in  practice,  especially  in  the  far  field  case  (Aki 
and  Richards,  1980,  p.123).  Also  note  that  the  desired  responses  for  P 
and  SV  excitations  could  be  obtained  by  an  appropriate  superposition  of 
the  responses  to  a  P-wave  source  and  to  a  mixed,  P-  and  SV-wave  source. 

Since  we  are  assuming  an  elastic  medium,  there  will  be  continuous 
conversion  between  P-type  and  SV-type  seismic  waves  as  the  inhomogeneous 
medium  is  penetrated.  This  makes  the  problem  far  more  complex  than  the 
acoustic  problem,  for  which  a  fast  algorithm  solution  has  already  been 
found  (Yagle  and  Levy,  1983).  There  may  also  be  propagation  of  SH 
waves,  if  the  impulsive  SV-wave  source  is  not  completely  polarized.  Since 
any  SH  waves  will  be  completely  decoupled  from  the  P  and  SV  waves,  these 
waves  will  not  be  considered  further  in  this  paper.  It  is  the  complexity 
of  wave  propagation  in  the  elastic  medium  that  allows  the  recovery  of  all  three 
medium  parameter  profiles  as  functions  of  depth  instead  of  travel  time. 

We  now  define  the  following  quantities: 


1/2 

a(z).  =  (XX  (z)  +  2]J  (z) )  /p  (z) )  =  local  P-wave  velocity  (la) 

B(z)  =  (y  (z)/p  (z) ) 1//2  =  local  S-wave  velocity  (lb) 

p  =  horizontal  ray  parameter  (1c) 

sin  0  (z)  =  a(z)p  =  sine  of  local  angle  between  P-wave 

ray  and  vertical  (Id) 

sin  0  (z)  =  $(z)p  =  sine  of  local  angle  between  S-wave 

ray  and  vertical  (le) 

a* (z)  =  a(z)/cos  0^(z)  =  local  vertical  P-wave  velocity  (If) 

3' (z)  =  3(z)/cos  0g(z)  =  local  vertical  S-wave  velocity.  (lg) 
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We  also  define  the  vector 


g  ( t ,  x ,  z ) 


u^ (t,x,z) 

u  (t,x,z) 
z 

Tzx(t,x,z) 

T  (t,X,z) 

.  zz 


(2) 


where  u  and  u  are  the  horizontal  and  vertical  components  of  the 
x  z 

displacement  and  where  x  and  x  are  the  horizontal  and  vertical 

zx  zz 

% 

tractions  on  an  element  perpendicular  to  the  z  axis. 

An  impulsive  plane  wave  bQ5 (t-px-qz)  is  used  to  probe  the  elastic 
medium.  Here  6(*)  denotes  the  Dirac  delta  function,  and  q  is  the  vertical 
ray  parameter  just  below  the  surface  (for  a  free  surface) ,  or  in  the 
homogeneous  half-space  above  the  medium.  The  Fourier  transform  of  this 
plane  wave  is  (bQ  exp- jtoqz)  exp-joopx.  Since  the  horizontal  ray  parameter  p 
is  independent  of  depth,  we  may  write  the  Fourier  transform  of  the  vector 
(2)  for  z>0  (inside  the  medium)  as 

A  A 

g(oo,x,z)  =  f  (w,z)exp-jcopx  .  (3) 

From  Aki  and  Richards  (.1980,  p.269),  the  propagation  of  seismic 
waves  in  an  inhomogeneous,  layered,  continuous  elastic  medium  is  described 
by 

8f/3z  =  A(z)f(to,z) 


(4) 
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where 


0 

-jcop 

i/y 

0 

--jcopA/  (A+2y) 

0 

0 

l/(A+2y) 

2  2  2 

4(0  p  y  (A+y)/(A+2y)-pco 

0 

0 

-  jtopA/  (A+2y) 

(5) 

0 

CN 

3 

CL 

1 

-j(0p 

0 

In  the  next  section  we  diagonalize  equation  (.4)  ,  defining  upgoing 
and  downgoing  P  and  SV  waves.  Appropriate  weightings  of  the  eigenvectors  of 
A(z)  will  be  necessary  to  put  equation  (4)  into  a  form  suitable  for  a 
fast  algorithm. 

TRANSFORMATION  OF  THE  PROPAGATION  EQUATION 
It  is  well-known  (e.g.  Claerbout  (1968) }  that  changing  variables 

A 

in  equation  (4)  from  f(co,z)  to  R(z)f((o,z),  where  R(z)  is  the  matrix  of 
row  eigenvectors  of  A(z) ,  diagonalizes  equation  (4)  into  upgoing  and 
downgoing  waves.  In  the  present  context  it  will  be  necessary  to  weight 
the  row  eigenvectors  of  A(z)  in  order  to  obtain  a  recursive  algorithm. 

Thus  we  define 

w((0,z)  =  X  (z)  R(z)  f  (0),  z)  (6) 

where  X  is  a  diagonal  matrix  whose  elements  weight  the  row  eigenvectors 
of  A(z) .  We  may  then  write 

f(.(0,z)  =  R  ^X  ^w((jO,z)  =  CX  ^w(W,z)  (7) 

where  C(z)  =  R(z)  1  is  the  matrix  of  column  eigenvectors  of  A(z) . 


Taking  the  partial  derivative  of  equation  (7)  with  respect  to  z  and 
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premultiplying  by  XR  yields 


3w/3z  =  [A-  (X(r3c/3z)X  1  +  X(3/3z(x  1)))]w 


(8) 


where 


A  =  RAC  =  diag  [- jw/a '  ,  —  jco/3  '  ju/a',  jw/g ' ] 


(9) 


We  now  choose  the  elements  of  the  diagonal  matrix  X  so  that  the  (diagonal) 
term  X3/3z(,X  )  =  -  (3/3z)log  |x|  zeroes  the  diagonal  elements  of 

X(r3c/3z).X  \  This  is  straightforward,  and  the  result  is 

X  =  diag  I  (ap  cos  0  )^2,  (gp  cos  0  (apcos:@  _)  ,  (gp  cos  0  )1//2]. 

P  s  id  s 

(10) 

We  recognize  the  components  of  X  as  the  square  roots  of  the  P-wave  and  SV-wave 

A 

impedances.  Hence  weighting  the  components  of  Rf  by  these  quantities  normalizes 
the  energy  fluxes  moving  upwards  and  downwards. 

Inserting  equation  (10)  in  equation  (8)  results  in 


-j0)/a'  -t 
t 


3w/3z 


-r 


-t 

c 

*-r 

P 

-r 

c 

-ju/8 ’ 

-r 

-r 

c 

s 

-r 

c 

jo)/af 

c 

-r 

s 

t 

c 

jco/3 

(ID 


where 


r  (z)  =  (l/2-2$2p2) (3/3z)log  p(z)  -  4g2p2(3/3z)  log  g(z) 


+  1/ (2-2a2p2) ( 3/3z)  log  a(z) 


(12a) 
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1/9  2  2  2 

rc(z)  =  -  (p/2)  (a'3’)  7  ((1-23  p  +  23  /ot'B')  (3/3z)log  p(z) 

-{462p2  -  432/a'B') (3/3z)  log  3(z) )  (12b) 

rg (z)  =  -(1/2  -  2g2p2) (9/3z)  log  p(z) 

-  (1/  (,2-2g2p2)  -  462p2)0/3z)  log  0(z)  (12c) 

tc(.z)  =  (p/2)  (a'3,)1/2((l-2B2p2  -  232/a'3')  (3/3z)  log  p(z) 

-(462p2  +  432/a'3") (3/3z)  log  3(z))  (12d) 


and  the  quantities  in  equations  (12)  have  the  following  interpretations: 


r  (z)  =  reflection  coefficient  for  a  reflected  P  wave 
P  generated  by  a  P  wave; 

r  (z)  =  reflection  coefficient  for  a  reflected  wave  generated 
by  a  wave  of  the  opposite  type; 

r  (s)  =  reflection  coefficient  for  a  reflected  SV  wave  generated 
by  an  SV  wave; 

t  (z)  =  transmission  coefficient  for  a  transmitted  wave  generated 
by  a  wave  of  the  opposite  type . 

We  have  used  here  notations  similar  to  those  of  Chapman  (1974)  and 
Kennett  and  Illingworth  (1981) .  The  physical  meaning  of  the  reflection 
coefficients  is  illustrated  in  Figure  1,  which  describes  an  infinitesimal 
section  of  a  lattice  filter  structure  which  implements  the  elastic  wave 
equation  (11)  .  Note  that  the  elementary  delay  elements  ==  exp  -  jcoA/a"  (z) 
and  Dg  =  exp  - j ojA/3 '  (z)  appearing  in  Figure  1  vary  with  depth.  The 
lattice  structure  of  Figure  1  can  be  viewed  as  a  generalization  of  lattice 
filters  used  in  speech  processing  (Markel  and  Gray,  1978)  and  linear 


estimation  theory  (Makhoul,  1977). 
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In  the  next  section  we  use  the  transformed  equation  (11)  to  obtain 
a  fast  inversion  algorithm. 


INVERSION  ALGORITHM 


Recall  that  the  first  experiment  consisted  of  probing  the  medium 

A 

with  a  planar  impulsive  P  wave.  Since  the  first  component  of  w(co,z) 
corresponds  to  a  downgoing  P  wave,  we  may  write  its  inverse  Fourier 
transform  w(t,z)  as 


w(,t,z)  = 


rb  6  (t-T  (z) ) 

P  P 

w1 (t,z) 

0 

+ 

w2 (t,z) 

0 

w3 (t,z) 

0 

w4(t,z) 

u  (t 


T  (z) 
p 


where 


T  (z)  =  /  d£/a'  (X) 

p  0 


denotes  the  vertical  travel  time  for  P  waves,  and 


u(t) 


1  for  t>0 
0  for  0 


(13) 


(14) 


(15) 


is  the  unit  step  function.  The  second  term  in  (13)  reflects  the  causality 
of  the  excitation:  There  can  be  no  wave  at  depth  z  until  the  excitation 
has  had  time  to  reach  depth  z. 

Taking  the  inverse  Fourier  transform  of  equation  (11) ,  inserting 
the  expression  (13) ,  and  equating  coefficients  of  6 (t-T  )  yields 

r  (z)  =  2w. (T  (z) ,z)/(a’ (z)b  )  (16a) 

p  3  p  p 

r  (z)  =  w.(x  (z), z)(l/a'(z)  +l/0'(z))/b  . 
c  4  p  p 


(16b) 
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Now,  for  the  second  experiment,  the  excitation  is  a  downgoing, 
impulsive  SV  wave.  Since  the  second  component  of  w(t,z)  corresponds  to 
such  a  wave,  we  have  for  this  experiment 


0 

W  (t,z) 

b  6(t-x  (z)) 

+ 

w„ (t,z) 

s  s 

2 

0 

w3 (t,z) 

0 

w4 (t,z) 

where  the  waves  w^(t,z)  have  the  form 


(17a) 


w^  It ,  z ) 


n. (t,z)u(t-T  (z) ) 
i  P 


+  q. (t,z)u(t-T  0 z ) ) 
X  s 


(17b) 


and  where  the  vertical  travel  time  for  SV  waves  has  been  defined  as 


T  (z)  =  .  ff  d£/B'(Ji)  .  (18) 

s  0 

Note  that  the  form  of  equation  (17)  differs  from  that  of  equation 
(13).  This  is  because  in  the  SV  experiment  the  impulsive  excitation  (an 
SV  wave)  does  not  coincide  with  the  wavefront  (a  P  wave) .  In  the  P  experiment 
both  of  these  were  P  waves  and  hence  coincided. 

Proceeding  as  above,  we  obtain  (q^ (z)  is  defined  in  equation  (17b)) 

rclz)  =  q3(Ts(z).,z)  (1/a'  (z)  +  l/0'(z))/bs  (19a) 

r  (z)  =  2q  (x  (z)  ,z)  /  (B’(z)b  )  .  (19b) 

s  4  s  s 

The  importance  of  equations  (16)  and  (19)  is  that  they  permit  computa¬ 
tion  of  the  reflection  coefficients  at  any  depth  provided  the  waves  w(t,z) 
are  known  at  that  depth. 
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Next,  equations  (12a-c)  are  written  as  a  matrix  equation: 


r  (z) 

P 

(3/3z) 

log  p(z)~ 

rc(z) 

=  M  (z ) 

(3/3z) 

log  8(z) 

_rs(Z)_ 

(3/3z) 

log  a(z) 

where 


(20a) 


M(z) 


2  2 

1/2  -  28  p 


„q2  2 

-48  p 


-£{.i-282p2  +  282/a,8')  £(482p2  -  4B2/a’B') 


-  d/2  -  2B2p2) 


-(1/(2-282P2)-482p2) 


with  £(z)  =  (p/2 )  (.a 1  8 1 ) i^2 .  Inverting  this  equation 


gives 


l/(2-2a2P2) 


0 

0 


(20b) 


(_3/3z).  log  p (z) 

r  (z) 
c 

(3/3z).  log  B  (z ) 

=  N(z)/m(z) 

r  (z) 
s 

(3/3z)  log  a(z) 

r  (z) 

— 

L  p  J 

(21a) 


where 


N(z)  = 


■(l/(2-282p2)  -  4B2p2)  -£(482p2  -  482/a'8')  0 


1/2  -  262p2 


_  q-a2p2)  (482p2-1) 


-£(1-282p2  +  282/a'8') 


2  2.  ,n2  2 


2(1-B2p2) 


and 


0 


-4Jl(l-a‘p‘)  (BY  +  82/a'8')  2m(l-a2P2) 

(21b) 


m(z)  =  (det  M(z) ) (2-2a2p2) 

=  £(1/2  -  3B2p2  -  B2/a'B'  +  2B4p4  +  2B4p2/a'B')/(l-B2p2)  .  (22) 
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Equations  (21).  function  as  update  equations  for  p(z),  B(z),  and  a(z). 

Note  that  p (z)  and  B(z)  are  updated  solely  from  r  (z)  and  r  (z) ,  and 

c  s 

then  r^(z)  is  used  to  update  a(z)  .  From  these  three  parameters,  any 
other  parameter  of  interest  (e.g.  A(z)  and  ]i(z))  may  be  quickly  found. 

Chapman  (1974,  p.  67)  gives  equations  similar  to  equations  (20);  however, 
Chapman's  equations  involve too  many  quantities  (y,p,8,a',  and  8")  and 
require  the  unobservable  transmission  coefficient  t  .  Thus,  they  are  un¬ 
suitable  as  update  equations. 

We  have  now  specified  all  of  the  equations  of  the  algorithm,  in 
differential  form.  The  algorithm  consists  of  equation  (11) ,  twice 
(one  for  the  experiment  involving  excitation  by  P  waves;  one  for  excitation 
by  SV  waves)  for  updating  the  up-  and  down-going  waves;  equations  (16)  and 
(19)  for  computing  the  reflection  coefficients;  equations  (21)  and  (22) 
for  updating  the  material  parameters  p(z),  B(z),  and  a(z)|  and  equation  (12d) 
for  computing  the  transmission  coefficient  tMz)  required  to  complete  the 
matrix  in  equation  (11) .  We  then  immediately  have,  for  each  z, 

y  (z)  =  82(z)p(z)  (23a) 

A(.z)  =  (a2(z)  -  2B2(z))p(z).  (23b) 

Next,  the  algorithm  is  discretized  in  order  to  clarify  the  recursions 
and  specify  in  what  order  quantities  should  be  computed. 

Discretization 

The  depth  coordinate  z  is  discretized  by  z=nA,  where  n  is  a  positive 
integer  and  A  is  the  discretization  length.  The  time  coordinate  t  is  similarly 
discretized  by  t  =  mAt,  where  At  is  the  discretization  time. 
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Initialization 

It  is  assumed  that  all  material  parameters  (X,  ]i  ,  p,  and  hence 
a,  0,  a',  and  0')  are  known  at  the  earth's  surface.  Since  we  are  assuming 
a  free  surface,  the  waves  at  the  surface  are  determined  by  measuring  the 
velocity  components  over  time,  for  both  the  P  and  SV  experiments. 

Recursion 

We  start  off  with  knowledge  of  a(z),  0  (z)  ,  p  (z) ,  a'(z),  B'(z)  as 
well  as  that  of  all  up-  and  down-going  waves  at  depth  z,  from  the  previous 
iteration.  Let  w^(t,z)  represent  the  waves  in  the  P-wave  source  experiment, 

g 

and  w  (t,z)  represent  the  waves  in  the  S-wave  source  experiment.  For 
convenience,  we  identify  the  dimensionless  quantities 

B(z)  =  02(z)p2  =  sin20s(z).  (24a) 

G(z)  =  02 (z)/a ' (z) 0 ' (z)  =  (l/2)sin  20  cot  0  ,  (24b) 

s  p 

Then,  taking  the  inverse  Fourier  transform  of  equation  (11)  and  employing 
a  simple  Euler-Cauchy  approximation  to  the  various  derivatives  in  the 
differential  form  of  the  algorithm  yields  the  following  recursive 
algorithm: 

1)  Computation  of  the  reflection  coefficients.  From  equations  (16)  and 
(19) , 

r  (z)  =  2w? (T  (z)  ,  z)/'(bA)  (25a) 

P  .3  p  ,.p 

r  (z)  =  2WP(T  (z),  z)/:'(bA)  (25b) 

c  4  P  p 

r  (z)  =  2(w®(t  (z)+,z)  -  w®(t  (z)“,z))/(b  A).  (25c) 

s  4  s  4  s  .23 

where  b  =  b  a’ (z)/A  and  b  =  b  01 (z)/A  are  the  strengths  of  the  discretized 
p  P  s  s 

continuous  impulses. 
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Upon  going  from  continuous  time  to  discrete  time,  the  continuous -time 

impulse  b_^<5(t)  becomes  a  discrete-time  impulse  of  hight  b^/D^ ,  where  is  the 

differential  delay  time  at  depth  z  for  wave  type  i  (see  Fig,  1) .  Since  the  impulse 

has  been  spread  out  over  the  time  interval  ,  its  height  must  be  b^/D^  in 

order  to  maintain  its  area  b.  .  For  a  P-+P  reflection  D  =  A/a  '  (z)  .  For  a  P->S 

1  p 

reflection  the  two-way  delay  is  D  +  D  ,  hence  the  one-way  delay  is  half  of  this, 

F  s 

or  (_A/2)  (1/a '  (z)  +  1/3'  (z)).  Equations  (16)  and  (19)  are  thus  modified  to  eq.  (25). 

2)  Computation  of  auxiliary  quantities.  From  equations  (22)  and  (24) , 


B  (z)  =  32  (z)  p2  (26) 

G (z)  =  32(z)/ci'  (z)B'  (z)  (27) 

m (z)  =  (1/2  -  3B  -  G  +  2B2  +  2BG)£/(1-B)  (28) 

t  (z)  =  -  £(1/2  -  3B  +  G  +  2B2  -  2BG) r  (z) / ( (l-B)m (z) ) 
c  c 

+  2Br  (z)/m(z)  (29) 

s 

where  £(z)  is  defined  as  above. 

3)  Update  of  material  parameters.  From  equations  (21) , 


p(z+A)  =  p(z)  -  p(z) ( (l/(2-2B)  -  4B)r  (z)  +  4£  (B-G)r  (z))A/m(z) 

c  s 

(30) 

3 (z+A)  =  3 (z)  -  3 (z) ( (2B-l/2)r  (z)  +  £(l-2B+2G)r  (z))A/m(z)  (31) 

c  s 

a(z+A)  =  a(z)  +  a (z) (1-a2 (z)p2) (2r  (z)  -  ( (2B-1/2) / (l-B)m(z) )r  (z) 

P  c 

-  £ (4 (B+G) /m(z))r(z))A  (32) 

2  2  1/2 

a' (z+A)  =  a (z+A) /(1-a  (z+A)P  )  7  (33) 

3' (z+A)  =  3(z+A)/(l-32(z+A)P2)1/2  .  (34) 

4)  Wave  update.  From  'the  inverse  Fourier  transforms  of  equation  (11) , 

w  (t+A/a'(z),  z+A)  =  w  (t,z)  -  (t  (z)w  (t,z)  +r  (z)w  (t,z) 

1  X  o 


+  r  (z)w  (t,z) ) A 
C  4 


(35a) 
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w  (t+A/8' (z) ,  z+A)  =  w  (t,z)  -  (-t  (z)w  (t,z)  +  r  (z)w  (t,z) 

C  J.  C  j 

+  r  (z)w.(t,z))A  (35b) 

s  4 

w_  (t-A/Ot'  (z)  ,  z+A)  =w„(t,z)  -  (r  (z)w.  (t,z)  +  r  (z)w  (t,z) 

3  3  p  1  c  2 

+  t  (z)w. (t,z))A  (35c) 

c  4 

w.  (t-A/B1  (z)  ,  z+A)  =  w.  (t,z)  -  (r  (z)w  (t,z)  +  r  (z)w  (t,z) 

4  4  c  i  s  2 

-  t  (z)  w_  (t ,  z) )  A  (35d) 

c  3 

and  these  same  recursions  are  used  for  both  w^(t,z)  and  wS(t,z). 

At  this  point,  we  have  obtained  p(z+A),  a  (z+A),  (3  (z+A),  a' (z+A) 

3' (z+A),  and  all  eight  waves  at  depth  z+A.  Hence  the  recursion  is 

complete.  Each  step  in  the  recursion  can  be  implemented  as  one  stage  or 

section  of  ;a ladder- type  filter,  which  can  be  regarded  as  a  more  complex 

version  of  the  lattice  filter  commonly  encountered  in  spectral  estimation 

theory.  A  typical  section  of  this  ladder  filter  is  illustrated  in  Figure 

1.  The  downgoing  P  and  SV  waves  at  depth  z  enter  the  filter  section 

at  the  upper  left,  interact  with  each  other,  are  reflected  (due  to  the 

inhomogeneity  of  the  medium) ,  and  exit  at  the  upper  right,  now  at  depth 

z+A.  Upgoing  P  and  SV  waves  undergo  a  similar  experience  in  the  lower 

half  of  the  filter.  Note  how  this  filter  illustrates  the  physical  meaning 

of  the  reflection  coefficients  r  (z) ,  r  (z)  and  r  (z) ,  and  of  the 

pc  s 

transmission  coefficient  t  (z)  . 

c 

The  recursions  of  the  waves  in  z  and  t,  given  by  equations  (35), 
are  slightly  complicated,  so  the  recursion  patterns  are  illustrated  in 
Figures  2a  and  2b.  We  start  off  knowing  the  waves  at  depth  z  for  all 
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time,  and  wish  to  find  the  waves  at  depth  z+A.  Although  the  simultaneous 
time  and  depth  updates  may  make  it  seem  as  though  information  at  early 
times  is  being  lost,  recall  that  by  causality  there  can  be  no  wave  at  depth 
z  until  the  initial  excitation  has  had  time  to  reach  depth  z.  Thus  there 
is  no  information  to  lose  at  the  early  times. 

The  algorithm  that  we  have  described  above  for  reconstructing  p(z), 

A(z)  and  y(z)  works  even  if  some  turning  points  exist  for  the  P  and  SV 
waves  propagating  through  the  elastic  medium.  However,  in  this  case 
p,  A  and  y  can  only  be  reconstructed  up  to  the  depth  z^  where  the  ray 
path  for  the  P  wave  becomes  horizontal.  Note  that  along  rays  associated 
with  the  P  and  SV  waves 

(36) 

so  that  unless  a(z)  <  1/p  for  all  z  (in  which  case  we  have  also  6(z)  <  1/p), 


sin  6  (z)/a(z)  =  sin  0  (z)/$(z)  =  p  =  constant 
P  s 


the  angle  0^(.z)  will  become  imaginary  at  some  depth  z  .  Physically,  this 
situation  results  in  evanescent  waves  where  the  waves  decay  exponentially 
with  depth.  This  causes  no  problem  in  the  reconstruction  algorithm  until 
z=z^„  at  which  point  a'  (z)  -*■  00 .  Then,  the  waves  w^(z,t)  and  wS(z,t) 
cannot  be  propagated  further,  and  the  material  parameters  are  reconstructed 
only  up  to  depth  z  . 

PHYSICAL  INTERPRETATION  OF  THE  ALGORITHM 


The  basic  principle  behind  the  algorithm  is  the  concept  of  layer- 
stripping,  which  is  discussed  in  Bruckstein,  Levy,  and  Kailath  (1983) . 
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At  each  depth  the  downgoing  and  upgoing  waves  are  being  scattered  (i.e., 
reflected  and  transmitted)  due  to  the  varying  material  parameters.  By 
including  an  impulse  in  the  initial  excitation,  the  reflection  coefficients 
may  be  measured,  since  the  impulse  is  an  easily  recognizable  tag.  The 
reflection  coefficients  are  then  used  to  update  the  material  parameters 
at  that  depth.  When  another  infinitesimal  layer  is  identified  it  is 
"stripped  away",  and  the  next  layer  is  examined  in  the  same  manner. 

The  layer- stripping  algorithm  is  closely  related  to  the  dynamic 
deconvolution  algorithm,  which  is  commonly  associated  with  inverse 
acoustic  problems  (.Robinson,  1982;  Yagle  and  Levy,  1983.;,  Bube  and 
Burridge,  1983)..  Both  algorithms  employ  up-  and  down-going  waves  being 
scattered  by  an  inhomogeneous  medium  with  the  material  parameters 
determined  from  reflection  coefficients.  However,  layer-stripping 
algorithms  may  take  a  wide  variety  of  forms,  while  dynamic  deconvolution 
is  generally  associated  with  a  specific  form  —  the  Schur  algorithm. 

The  waves  (.elements  of  w).  used  in  the  algorithm  are,  in  the  Fourier 
transform  domain 

T  /Z  1/2  +  jWZ  1/2U  (37a) 

P  P  P  P 

T  /Z  1/2  +  jtoZ  1/2U  (37b) 

s  s  s  s 

where  the  upper  sign  is  used  for  upgoing  waves  and  the  lower  sign  for 

downgoing  waves,  the  impedances  Z  and  Z  are  defined  by 

P  s 

Z  (z)  =  a(z)p(z)  cos  0  (z)  (38a) 

p  p 

Z  (z)  =  $(z)p(z)  cos  6  (z)  ,  (38b) 

s  s 

AAA  A 

and  T  ,  T  /  U  and  U  are  the  stresses  and  displacements  along  the  ray  path 

jp  Q  jp  CJ 

for  P  waves,  and  perpendicular  to  the  ray  path  for  SV  waves. 
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This  shows  that  the  up-  and  down-going  waves  have  been  normalized 
at  each  depth,  so  that  the  energy  in  a  wave  is  the  square  of  its  amplitude, 
regardless  of  the  medium  around  it.  This  preserves  the  equality  of 
incoming  and  outgoing  energy  fluxes  in  any  region.  Note  that  in  a  homogeneous 
medium  the  waves  (37).  become  simply  the  energy-normalized  velocities,  which 
were  the  state  variables  used  by  Shiva  and  Mendel  (1983)  for  a  discrete 
(each  layer  is  homogeneous)  layered  elastic  medium. 

It  should  also  be  noted  that  if  the  medium  is  discretized,  i.e. 
modelled  as  a  welded  stack  of  thin,  homogeneous  layers  with  material 
parameters  varying  only  between  different  layers,  then  r3c/3z  may  be 
interpreted  as  a  scattering  matrix  for  the  layer  at  depth  z .  To  see 
this,  replace  (3/3z)  log  p(z)  =  (3/3z) p (z) /p (z)  by  the  discrete  approxi¬ 
mation  Ap(z)/p(z),  and  do  the  same  for  B(z)  and  a(z)  .  Then  equations 
(.12)  become  the  reflection  and  transmission  coefficients  at  an 
interface  (Aki  and  Richards,  1980,  p.  153).  Thus  discretization  of  the 
algorithm  is  equivalent  to  a  physical  discretization  of  the  medium. 

RESULTS  OF  A  COMPUTER  RUN  OF  THE  ALGORITHM 

The  algorithm  was  tested  by  running  it  on  the  synthesized 
impulse  response  of  a  twenty-layer  medium.  The  variation  of  medium 
parameters  from  one  layer  to  another  was  made  small  (around  2%) ,  in 
order  to  simulate  a  continuous  layered  medium.  This  is  important,  since 
the  differential  updates  assume  a  continuously  varying  medium;  the  algorithm 
cannot  handle  sharp  changes  in  medium  properties  unless  the  step  size  A  is 
made  smaller  in  such  regions.  The  medium  velocities  and  step  size  A  were 
scaled  down  by  a  factor  of  1000,  so  A  =  0.1m  instead  of  100m. 
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The  response  of  the  medium  to  impulsive  plane  P  and  SV  waves  was  generated 

in  the  frequency  domain  using  the  reflectivity  method  (Aki  and  Richards,  1980, 

p.  393) .  A  FORTRAN  program  given  by  Kind  (1976)  was  used  to  compute  the  plane 

wave  transfer  functions  R  ,  R  ,  and  R  at  512  frequency  points  (integer 

pp  ps  ss 

multiples  of  0.78  Hz).  Each  of  these  was  divided  by  j2irf  and  a  discrete  inverse 
Fourier  transform  taken.  This  synthesized  sample  Step  responses;  taking  diffe¬ 
rences  and  dividing  by  the  discretization  time  At  =  0.0053  yielded  the 
discretized  impulse  responses.  It  should  be  noted  that  careful  attention  must 
be  paid  to  signs  in  going  from  potential  reflection  responses  to  velocity 
reflection  responses;  see  (Aki  and  Richards,  1980,  p.  191). 

The  impulse  responses,  scaled  by  1/At  for  convenience,  are  plotted  in 
Figure  3.  Although  the  responses  were  computed  for  t  =  0  up  to  t  =  2.565 
to  avoid  aliasing  problems,  the  responses  beyond  t  =  1.3s  were  essentially 
zero  and  are  not  shown.  Note  that  the  peaks  corresponding  to  strong  primary 
reflections  are  smeared  out.  This  is  due  in  part  to  the  use  of  a  DFT ,  which 
in  this  case  is  tantamount  to  bandpass-filtering  the  data  with  a  filter  with 
pass  band  0.78  Hz  -  400  Hz,  Since  the  strengths  of  the  primary  reflections  are 
especially  important  to  the  algorithm, this  smearing  might  be  expected  to  hamper 
its  performance.  However,  this  evidently  did  not  happen. 

The  impulsive  plane  wave  responses  were  then  used  to  initialize  the 
upgoing  P  and  SV  waves,  and  the  algorithm  was  run  on  a  VAX7II/782  computer. 
Results  are  shown  in  Figure  A,  It  can  be  seen  that  the  agreement  between 
the  actual  and  algorithm-generated  medium  parameter  profiles  is  quite  good, 
with  less  than  5%  error  everywhere. 

It  should  be  noted  that  the  algorithm  was  not  tested  under  perfect  con¬ 
ditions.  Bandlimiting  of  the  frequency  response  resulted  in  the  time  response 
being  smeared  over  two  or  three  samples,  and  the  medium  itself  was  discrete, 
so  that  some  error  may  be  expected  in  the  update  equations.  Nevertheless, 
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the  algorithm  performed  quite  well.  Small  amounts  of  additive  noise  in 
the  data  may  be  handled  by  running  the  algorithm  concurrently  at  several 
different  values  of  the  stacking  parameter  p,  averaging  the  updated  medium 
parameters  at  each  depth,  and  using  these  averaged  updates  in  the  algorithms. 
Sharper  variations  in  the  parameter  profiles  could  be  handled  by  temporarily 
reducing  the  step  size  A  if  the  reflection  coefficients  get  too  large;  this 
would  also  reveal  the  sharp  variation  in  more  detail.  These  possibilities 
are  subjects  of  current  research. 

CONCLUSION 

A  fast  algorithm  for  recovery  of  material  parameter  profiles  as 
functions  of  depth  for  the  case  of  seismic  wave  propagation  through 
a  continuous  elastic  medium  has  been  given.  The  algorithm  has  been  specified 
both  in  differential  form  and  in  a  simple  discretized  version  that  details 
its  recursive  nature.  The  algorithm  works  on  a  layer-stripping  principle, 
and  appears  to  be  much  faster  and  easier  computationally  than  previous 
solutions  to  the  inverse  seismic  problem  for  an  elastic  medium.  A  physical 
interpretation  of  the  algorithm  was  also  discussed  in  terms  of  a  lattice 
filter  concept  showing  how  the  first  reflection  of  various  wave  types  at  each 
depth  yields  the  medium  reflection  coefficients  at  that  depth,  from  which 
the  medium  parameters  at  that  depth  can  be  differentially  updated  to  the 
next  depth.  Results  of  a  computer  run  of  the  algorithm  on  the  impulsive 
plane  wave  response  of  a  twenty-layer  medium  show  that  the  algorithm  works 
satisfactorily. 
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More  work  needs  to  be  done  in  investigating  the  speed,  stability, 
and  performance  of  this  algorithm  on  real-life  data.  The  effects  of 
noise  and  modelling  errors  on  this  algorithm  should  also  be  the  subject 
of  further  research. 
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FIG.  2b 


Recursion  pattern  for  updating  the  upgoing  waves 


Fig.  3a.  P->P  impulse  response,  scaled  by  1/At  for  conven 


P-KP  impulse  response/At 


Pig.  4a  Comparison  between  actual  and  computed  P  wave  speed  profiles 
(2  =  actual,  3  =  computed) .  Both  profiles  and  depth  have  be 
scaled  by  1000. 
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